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We study the relaxation and damping dynamics of an ultracold, but not quantum degenerate, 
gas consisting of dipolar particles. These simulations are performed using a direct simulation Monte 
Carlo method and employing the highly anisotropic differential cross section of dipoles in the Wigner 
threshold regime. We find that both cross-dimensional relaxation and damping of breathing modes 
occur at rates that are strongly dependent on the orientation of the dipole moments relative to the 
trap axis. The relaxation simulations are in excellent agreement with recent experimental results 
in erbium. The results direct our interest toward a less-explored regime in dipolar gases where 
interactions are dominated by collision processes rather than mean-field interactions. 

PACS numbers: 


I. INTRODUCTION 

Much of the attention on ultracold dipolar gases has 
heretofore focused on the quantum degenerate regime, 
where dipolar interactions can significantly influence the 
behavior of the gas through the mean-field. Aspects of 
this influence include changing the shape and mechanical 
stability of the gas M , as well as altering the excitation 
spectrum to include low-energy roton modes in a Bose- 
Einstein condensate M. A host of related phenomena 
have been predicted and observed [l0l4l3| . driven by the 
direct action of the long-ranged, anisotropic dipolar in¬ 
teraction on the particles’ motion. 

By contrast, gases at a slightly higher temperature 
behave more classically, and their mean-field energy is 
overcome by kinetic energy as the prime source of dy¬ 
namics in the gas. In such a situation the strength and 
anisotropy of the dipolar interactions can be made man¬ 
ifest through collisions, rather than through mean-field 
effects [l4j . A very recent experiment showed this ex¬ 
plicitly, finding that collisional relaxation of a gas of er¬ 
bium atoms at ~ 400 nK occurred on time scales that 
varied by a factor of four, depending on the orientation 
of the atoms’ magnetic dipole moments 0 ]. This land¬ 
mark result illustrates the potential for anisotropic dipo¬ 
lar scattering to profoundly influence the kinetics of a 
cold, thermal gas, from rethermalization and relaxation, 
to viscosity and the propagation of sound, to name but 
a few features. 

In this article we construct a model of the cold, non¬ 
degenerate dipolar gas by numerically solving the Boltz¬ 
mann equation. The model is based on the direct sim¬ 
ulation Monte Carlo (DSMC) algorithm 0,111 , which 
is appropriate to the dilute limit found in experiments, 
when the mean-free path A m f of the atoms in the gas is 
comparable to, or larger than the characteristic scale L of 
the gas (i.e., Knudsen number Mi = A m f /L >1). Using 
this model, we explore the thermal relaxation and damp¬ 
ing of a dipolar gas that is suddenly taken out of equi¬ 
librium. Where applicable, our results are in excellent 
agreement with the return to equilibrium of the erbium 
gas in Ref. jl5| , and in particular describe the dependence 


of relaxation rate on polarization direction of the dipoles. 
Further, we characterize the damping rate of breathing 
mode oscillations generated in the gas, finding that this 
damping is also strongly dependent on polarization, and 
is slower than the rethermalization rate. We also eval¬ 
uate the relevance of mean-field interactions in the gas. 
Although the density of erbium in the cross-dimensional 
relaxation experiment was not sufficiently high to observe 
mean-field effects, we briefly discuss how to modify the 
DSMC method to include such physics (using particle- 
in-cell methods for a dipolar-Vlasov equation). 

The outline of this paper is as follows: In section HT1 we 
introduce and discuss the details of a cross-dimensional 
rethermalization experiment which we will model. In sec¬ 
tion [ITT] we provide a very brief introduction to the Boltz¬ 
mann equation and discuss its historical significance in 
statistical mechanics. Section UlI Bl outlines the basic fea¬ 
tures of our DSMC algorithm, and section [TiTC] discusses 
the differential scattering cross sections for low-energy 
dipolar interactions. Section m discusses and quantifies 
the mean-field interaction in the gas. Section fVl contains 
our results for fermions, and compares these results to 
experimental data. Section eh contains similar results, 
but for bosons. In section IVIII we draw conclusions and 
discuss possible avenues for future research. 


II. CROSS-DIMENSIONAL RELAXATION OF A 
DIPOLAR GAS 

For concreteness, we here contemplate the experimen¬ 
tal situation of Ref. [TB] . We employ the notation of that 
experiment, and use the same values of trap frequencies, 
density, and species (erbium). We stress, however, that 
the simulations can be made completely general for cold 
dipolar gas experiments in the thermal regime, including 
polar molecules. 

Experiments involving cross-dimensional relaxation 
have a long history in cold atoms, going back to the work 
with caesium! 18|. Other experiments include work on 
Bose-Fermi 0 ] and Fermi-Fermi mixtures 0 ]. The ex¬ 
perimental scenario we consider is shown in Fig. [T] The 
gas begins in the equilibrium state of an approximately 
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cylindrically symmetric trap, with the dipole alignment 
direction in the y—z plane of the laboratory reference 
frame. The gas is weakly trapped in the y direction, 
and tightly trapped in the x and z directions. The dipole 
alignment direction, e, makes an angle /3 with the y— axis. 

Over a (fast) time scale t ramp , the trapping frequency 
along the y —axis is significantly increased, sending the 
system out of equilibrium. The atoms, whose distribu¬ 
tion is initially still elongated along the y direction, gain 
extra momentum along this direction (over the time-scale 
of a quarter trap-period). Rethermalisation requires the 
redistribution of this additional momentum/potential- 
energy in the y direction into the x and z directions. 
Due to the highly anisotropic nature of the dipole-dipole 
interaction, the rate at which this rethermalization (re¬ 
distribution) occurs depends strongly on the angle, /3, 
between the dipole alignment direction and the y —axis 
(see Fig.[T]). 

This experiment was recently performed in Inns¬ 
bruck [ij| , as a very beautiful demonstration of the stan¬ 
dards in precision and control over cold-atomic systems. 
The atomic species used was 167 Er (a fermion), which 
has an exceptionally large magnetic dipole moment of 
7/j.b where yB is the Bohr magneton (compared to s 'Rb 
with 1/is and 52 Cr with 6/zb, 164 Dy has 10/zb). The ex¬ 
periment began with an initial temperature of 426nK. 
Relative to the density of the system, this corresponds 
to a regime; n\j* ss 0.25, where n is the average density 
in the trap, and At = hy/2 -k/{ rnksT) is the thermal de 
Broglie wavelength. In this sense, the system (although 
cold) is not deeply within a regime of quantum degen¬ 
eracy. This then implies that the classical Boltzmann 
equation should provide the appropriate theoretical de¬ 
scription. This being said, quantum-mechanical effects 
may indeed be a source of error in our simulations, and 
we attempt to quantify this in Section [V Dl 

In spite of this (relatively) low phase-space density, 
the system is still sufficiently cold such that the ratio 
between the thermal de Broglie wavelength, and a char¬ 
acteristic dipole-length scale; ad = Cdd’m/(8'rrh 2 ), (where 
Cdd = yoy 2 for magnetic dipoles, and Cdd = d 2 /e o for 
electric dipoles, yo and eg are respectively the permeabil¬ 
ity and permittivity of the vacuum, fi and d are the mag¬ 
netic and electric dipole moments) is; Ax/cid ~ 39. From 
this, we conclude that the two-body scattering physics is 
strongly within the quantum regime, and the differential 
scattering cross-sections are chosen accordingly fl4 ]. 


III. THE BOLTZMANN EQUATION FOR 
DIPOLAR GASES AND THE DSMC METHOD 



FIG. 1: (Color online) The initial state of the gas is shown 
in (a), where the atoms occupy the equilibrium state of an 
approximately cylindrically symmetric trap, elongated in the 
y direction. The dipole alignment direction is given by e. 
As depicted in (b), the experiment begins when the trap- 
frequency in the y direction is suddenly increased, sending 
the system out of equilibrium. Rethermalisation dynamics 
depends on the angle /3 between the dipole alignment direction 
and the y axis. 


A. General considerations 

The ability to trap and cool atoms with large mag¬ 
netic dipole moments, such as chromium [2l|, (22|, dys¬ 
prosium [23j, Q , and erbium [25] provide exciting pos¬ 
sibilities for observing novel many-body states (for a re¬ 
cent example, see Ref. [26}). Dipolar molecules are an¬ 
other source of_ potentially even stronger interactions in 
dipolar gases [13} l27H3l} . Developing theoretical tools to 
understand dipol ar g ases is currently a very active area 
of research ist tl3 [37f . 

A particularly challenging task in many-body physics 
is to develop theoretical methods for treating out-of¬ 
equilibrium physics. To this end, we report on our 
progress towards a general tool for simulating out-of- 
equilibrium dynamics of the normal dipolar gas. Our 
approach is based on the Boltzmann equation, which we 
solve using the DSMC algorithm. The motivation for us¬ 
ing DSMC typically occurs when the assumptions of fluid 
mechanics (which generally centres around some form of 
the Navier-Stokes equation) break down, and one must 
account for the granular nature of matter (usually, al¬ 
though not exclusively, via statistical mechanics). Bird’s 
DSMC algorithm has evolved over recent decades into a 
remarkably versatile and useful tool which has been ap¬ 
plied across seemingly disparate fields of research [3SJ, [39} . 

Stochastic particle methods, such as DSMC, have been 
applied to ultra-cold gases in a number of previous works. 
For instance, a variation of the method we describe here 
was used to study evaporative cooling enroute to Bose- 
condensation |40j. In Ref. [4l|, collisions between two 
thermal clouds near a d -wave resonance was simulated. 
The results compared very nicely to experiment |42} . 
Other examples include the s tudy of collective modes in 
finite temperature dynamics [431 - 145} . sympathetic cool¬ 
ing of molecules [46}, and degenerate Fermi gas dynam¬ 
ics [47j . To our knowledge, our work is the first time 
dipolar differential cross sections have been used 0- 
This reduces the efficiency of the DSMC by introducing 
a rejection-sampling algorithm to sample the differential 
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cross sections. However, we find that numerical conver¬ 
gence is still easily attainable on standard commodity 
hardware. 

The classical Boltzmann equation describes the statis¬ 
tical mechanics of particles in a many-body system with 
two-body elastic collisions. Its modern derivation typi¬ 
cally involves truncation of the BBGKY hierarchy [4|| 
such that two-body (and higher) distribution functions 
factorize into products of single-body distribution func¬ 
tions (this assumption was referred to by Boltzmann 
as the stosszahlansatz , or the assumption of molecular 
chaos). The equation for a single component gas reads 


d_ 

dt 


^•V r 

m 


■F.Vr 


/ = C[f] 


(1) 


where / = /( r, p; t) is the single particle phase-space dis¬ 
tribution, i.e. fd 3 rd 3 p is the expected number of atoms 
within the phase-space volume (r, p) —y (r+d 3 r, p+d 3 p), 
m is the particle mass, F denotes the external forces act¬ 
ing on the system, i.e F = — V r C/(r,t) where U(r,t) is 
some external potential (trapping potential), and finally 


c|/1 = /^/‘ in ^ lp - p ‘ ll/ ' /: - // ' 1 (2) 
is the collision integral. We have used the common no¬ 
tation = /(r, p^;t). In principle, one may wish to 
include a mean-field contribution into the external poten¬ 
tial. We discuss the relative importance of this mean-field 
term in section HVl and demonstrate its insignificance for 
the purpose of simulating the experiment in Ref. fl5j . 

The collision integral in Eq. 0 provides a mechanism 
for rethermalization via two-body collisions. Two parti¬ 
cles (coinciding at the point r) collide with momenta p 
and pi, and emerge from the collision with momenta p' 
and p^. Net energy and momentum are conserved in the 
collision, meaning 

P = P' (3a) 

I Prel | = | Prel I ( 3b ) 

where = (p (b + p^)/2 and p^ = 

denote center-of-mass, and relative, momentum respec¬ 
tively. The differential cross-section 

srK< p-p-) < 4 > 

contains information regarding the likelihood of two par¬ 
ticles colliding (given an incident relative momentum 
p re i), and the likelihood of a post-collision relative mo¬ 
mentum given by p( el . Intriguingly, cross-sections which 
exhibit time-reversal symmetry, that is 4^(p re i, p( el ) = 
jifi (p( e i■ Prel) 1 yield irreversible dynamics in the Boltz¬ 
mann equation as demonstrated by Boltzmann’s famous 
H-theorem [4!|. The relevant differential cross-section 
for dipolar particles has been derived and discussed in 
detail in a recent article [l4| , and we will briefly sum¬ 
marise the necessary results in section Ull Cl 


Analytic solutions to the Boltzmann equation are dif¬ 
ficult to come by [49|, 0]. An important exception 
are the well known equilibrium solutions, the Maxwell- 
Boltzmann distribution; 


N 

/(r, p; t) = /mb(i% P) := ^ exp 


p 2 /2 m + U{ r) 


knT 


(5) 

where ks is Boltzmann’s constant, T is the temper¬ 
ature, N is the total number of particles, and Z = 


J d 3 rd 3 p exp 


p 2 /2m+U(r) 

k B T 


gives the correct normalisa¬ 


tion. Using the conservation laws in Eq. ©, it is straight¬ 
forward to see that C[/mb] = 0, and /mb is a stationary 
solution to Eq. ©■ 

We wish to solve the Boltzmann equation © under the 
following dynamical scenario: Starting from an equilib¬ 
rium initial distribution, Eq. ©, we change the trapping 
potential U( r, t) = \m [co 2 x 2 + uj y (t) 2 y 2 + w 2 z 2 ], where 


Uy(t) = < 


,(0) 


■ UJ- 


( 0 ) 


, \/l + ■ 




t < 0 

0 t 1 t ram p 

t . ' t ram p 


( 6 ) 


such that, over the ramp-time f ramp , the trap frequency in 
the y-direction is changed by a factor y/1 + s. The choice 
of square-root dependence on time corresponds to lin¬ 
early increasing the laser power in an optical dipole trap. 
The spatial anisotropy, created by the dipole-alignment 
direction, creates a bias for scattering into particular mo¬ 
mentum states. This has implications for the rate of 
rethermalization, which becomes dependent on the angle 
between the y-axis and the dipole alignment direction. 
We investigate the rethermalization dynamics as a func¬ 
tion of this angle, and compare it to the experimental 
work of Ref. [Tf|. 

Attempting to solve the Boltzmann equation by dis¬ 
cretizing the temporal axis and the phase-space dimen¬ 
sions is a futile exercise as (for all but the most triv¬ 
ial cases) one will always run out of computational re¬ 
sources, before numerical convergence is achieved. Vi¬ 
able alternatives to find an approximate solution in a 
close-to-equilibrium scenario do exist however. For in¬ 
stance, the so-called method of moments approach was 
used in Ref. .36j to study collective excitations of two- 
dimensional dipolar fermions in a perturbative limit. In 
Ref. git a variational method was employed to predict 
relaxation behaviour in s-wave interacting gases. Our 
method is more generally applicable to a wider variety 
of far-from-equilibrium scenarios, although it is more nu¬ 
merically intense than other methods. 


B. The DSMC method 

The starting point for DSMC approximates the distri¬ 
bution function /, by N t test-particles each with posi¬ 
tion and momenta [iy,pi] which are found by randomly 
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sampling /(r, p; t = 0). That is 
n t 

/(r,p;0)«^^^(r-r i )5(p-p i ) (7) 

i=l 

where £ = N/N x is the ratio of real-particles to test- 
particles. The goal is to force the test-particles to evolve 
in time [rj(i), p;(t)] such that their relationship to / 
shown in Eq. 0, remains true at all times. The compu¬ 
tational complexity thus increases with Nt- 

On time scales, At, much shorter than the mean- 
collision-time, the evolution of each test particle is given 
by its classical trajectory in the potential. Assuming At 
is also much shorter than the trap period, this is well 
approximated using a predictor-corrector (symplectic in¬ 
tegrator) method, 

q* = Ti(t) + ^-Pi(i) (8a) 

2m 

Pi{t + At) = p*(t) + F z At (8b) 

rj(t + At) = qj + -—Pi(t + At), (8c) 

where Fj = — V qj 17(qi,t) is the external force acting on 
the i- th test-particle. This is often referred to as the free- 
streaming dynamics. Note that, if the classical trajectory 
of a single particle in the trap can be solved analytically 
(which is obviously straight forward in the case of a har¬ 
monic potential), then Eqs. (JHJ) can be replaced by this 
analytic solution. This provides an advantage in that At 
need not be small compared to the trap period (but still 
must remain small compared to the mean-collision-time). 
In effect, Eqs. © account for the left hand side of the 
Boltzmann equation as shown in Eq. ©■ 

In order to include the effects of the collision integral 
[on the right hand side of Eq. ©], a spatial grid is intro¬ 
duced, and the test-particles are binned into the volume- 
elements AV of this grid. This grid needs to be cho¬ 
sen carefully. The size of the volume-element effectively 
represents the finite resolution of the delta-function in 
our numerics. For this reason it needs to be small since 
all physical quantities will be coarse-grained over this 
volume-element. However, we will use the population 
of test-particles within each volume-element to stochas¬ 
tically check for collisions, and therefore, the volume- 
element must be large enough to contain multiple test- 
particles (in order to obtain reliable statistics). Being 
certain that one has the necessary combination of large- 
enough Nt and small-enough AV is an important nu¬ 
merical convergence test. 

Once the spatial grid has been established, we check 
N V (N V — l)/2 pairs of particles within the zzth volume 
element (7V„ is the population of the zzth volume ele¬ 


ment). In this step, the computational complexity ac¬ 
quires a N1 dependence, and simulations will become un¬ 
feasible if individual volume elements contain too many 
test-particles. The collision probability is given by 

Pij = C-^ylPrelKPrel) (9) 

where p re i = p, (t) — p, (tjiand 

/ dcr 

d0p -‘d^( Prel,P " el ) ( 10 ) 

is the total cross section (as a function of relative momen¬ 
tum betwen particles i and j), found by integrating the 
differential cross section over all solid angles of scattered 
relative momentum. Computational parameters must 
be chosen such that Pij -C 1. The collision proceeds 
if R < Pij where R is a randomly generated number, 
with uniform distribution between 0 and 1. If the colli¬ 
sion proceeds, we establish the post-collision relative mo¬ 
mentum p[. el by treating the differential cross-section ^ 
as a probability distribution for p[. el , and stochastically 
sample it using a rejection-sampling algorithm (see Ap¬ 
pendix [A] for more details). The center-of-mass momen¬ 
tum is conserved during the collision. In this way, at each 
time-step in our simulation, collisions are stochastically 
implemented, in correct accordance with the total-cross 
section, the local density, the local velocity distribution, 
and the differential scattering. 

Our numerical algorithm described here, has some sub¬ 
tle inferiorities compared to certain other algorithms de¬ 
scribed in the literature. Deficiencies include the absence 
of locally-adaptive spatial grids (to efficiently account for 
dramatic variations in spatial density), scaled collision 
probabilities (without which the number of operations in 
the algorithm scales as ~ N^, rather than a potential 
~ Nt scaling), and locally adaptive time steps [4l|, [52| . 
However, the cold atomic vapours under current consid¬ 
eration have relatively small numbers of particles, and 
we have thoroughly tested for, and found, excellent nu¬ 
merical convergence in all of our simulations. For this 
reason, we do not implement the complete set of modern 
sophistications within the DSMC. 


C. Differential scattering in dipolar gases 

The cross-section formulae used in this work were de¬ 
rived in Ref. (lij using the Born approximation for the 
scattering amplitude between two dipolar particles, with 
dipole moments aligned along an alignment direction e 
(we use ~ to denote a unit vector). The formulae are 

^^(Prel, p'el) = ISF,B(Prel, p( el )| 2 , where 
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3F(Prel,Prel) 
1?B(Prel; P r el) 


1 4(p rel .e)(p' el .£) - 2 [(p r cl-e) 2 + (p/l^) 2 ] (Prel-Prel) 


1 

a/2 L " a d 


1 - (Prel-P'el) 2 


o a 2 (Prel-e) 2 + 2 (p; el .g) 2 - 4(p re l.£)(p; el .£)(p rel .p; el ) , 4 


1 ^ (Prel.p/j) 2 


( 11 a) 

(lib) 


ad is the dipole length scale given by ad = ro/xo/i 2 / (871A 2 ), 
/io is the vacuum permeability, fi is the atoms magnetic 
dipole moment (n = 7//b in the case of erbium), and a is 
the s-wave scattering length. The subscripts F and B re¬ 
spectively correspond to fermionic and bosonic symmetry 
constraints ( 167 Er which was used in the experiment [15[ 
is fermionic). 

The total cross section, which we use to evaluate the 
collision probability in Eq. © can also be evaluated an¬ 
alytically 0, 

cr F (Prei) = [ 3 + 18cos 2 (?7) - 13cos 4 (?7)] (12a) 

CTB(Prei) = a 2 ^{72a 2 - 24a [l - 3cos 2 (ry)] 

+ 11 — 30 cos 2 ( 77 ) + 27 cos 4 ( 77 ) j (12b) 

where 77 = cos _1 (p re i.£) is the angle between the dipole 
alignment direction and the incoming relative momen¬ 
tum. Equation (fl2l) [(a) or (b) depending on whether 
the collision pair are identical fermions or bosons] is used 
in Eq. © to evaluate the likelihood of a collision. 

Once it has been established whether or not the col¬ 
lision occurs, the post-collision relative velocity is found 
by sampling the distribution function 


To convert between the lab-reference-frame and the 
collision-reference-frame, we find 

ef = [cos(7) cos(/ re i) cos(0 re i) - sii/y) sin(/ rel )] e 1 / 

+ [cos(7) sin(/ re i) cos(0 re i) - sin(7) cos(/ re i)] elf 
- cos(7) sin(0 re i)e3 ( 14 ) 

e 2 f = [- sin(7) cos(/ re i) cos(6>rei) - cos(7) sin(/ re i)] e“ 

+ [cos(7) cos(/ re i) - sin(7) sin(/ rel ) cos(0 re i)] eif 
+ sin(7) sin(0 re i)e3 ( 15 ) 

ef — cos(/ re i) sin(0 re i)e 1 / + sin(/ re i) sin(0 re i)e2-|- 

cos(6> rel )e3 ( 16 ) 

where the angle 

7 =acot{cos(0 re i) cot(/ E - cj) ie 1) 

- cot( 0 E )csc(</> E - /rei) sin( 0 re i)} ( 17 ) 

and 

Prel = sin( 0 re i) cos^re/e 1 / + sin( 0 re i) sin(/ rel )e2 

+ COs(6»rel)e3 , ( 18 ) 

£ =sin(0 E ) cos(/ E )e 1 1 f + sin(0 E ) sin(/ E )e2 + cos(6 | E )e3 . 

( 19 ) 


-PF,B(0rel, 4>rel', V) 


1 

°F,B(Prel) dtl 


(Prel, Prel) sill ^rel* 


(13) 

Note that we only need to sample 0 re i and / re i since 77 
is given to us by the (already known) incoming relative 
momentum of the collision pair. The collision-reference- 
frame (x c f, y c [, z c f) is defined such that the z c f-axis points 
along the direction of p rc i, and the dipole-alignment di¬ 
rection £ lies in the x c f z c { plane. The purpose of defin¬ 
ing, and operating within the collision-reference-frame is 
to make the analytic formulae of Eqs. (Illall and (lllbll as 
wieldy as possible. The coordinates 6 and 4> in Eq. m 
are the polar and azimuthal angles (respectively) of p( el 
in the collision-reference-frame. We (arbitrarily) decide 
to include the factor sin 6 into the definition of the proba¬ 
bility distribution function (rather than the metric) such 
that J^dcj) dd P f ,b( 0, <(>; ??) = 1- Sampling the prob¬ 
ability distribution in Eq. m is not simple, so we use 
a rejection sampling algorithm which we describe in Ap¬ 
pendix [A] 


We have used the common notation where e 1 /’^ 3 denote 
the standard (unit) basis vectors of Euclidean space in 
either the lab- (If) or collision- (cf) frame. The symbols 
6 e and </> E refer to the azimuthal and polar angles (re¬ 
spectively) of the dipole alignment direction in the lab 
frame, as shown in Eq. m- 


IV. DISCUSSION ON THE MEAN-FIELD 
INTERACTION 


In a more general situation the inclusion of a mean- 
field interaction may be desirable (HU, [54| . This requires 
an alteration to the Boltzmann equation m such that 
F = -Vl/(r, f) now consists of two parts, U(r,t) = 
f/ ex t(r,t) + U mf (r,t), an external potential U ex t and 
a mean-field potential U m {. Such an approach may 
be dubbed a dipolar-Vlasov equation in recognition of 
its similarity to the Vlasov equation used in plasma 
physics [55j. The mean-field potential is a dynamical 
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variable (away from equilibrium) found from the convo¬ 
lution 

Umf{r,t) = j d 3 r'n(r',t)V dd (r-r') (20) 

where n( r, t) = f d 3 p /( r, p; t) is the spatial number den¬ 
sity and Vdd(r) is the dipolar interaction between two 
particles separated by r. This is given by 


are difficult in this case. Although a numerical solution is 
not difficult, it only changes the result by a factor of order 
unity, and is therefore not of interest to us at this stage. 
The wonderfully elegant Fourier transform of Vdd(r) al¬ 
lows for the analytic calculation of e m f 0 


- 4 = 

48Vtt 3 \ a z ) 


(24) 


Vdd(r) 


Cdd 1 ~ 3(r • ef 

47 r r 3 


( 21 ) 


In general it is certainly true that the physics associated 
with the mean-field interaction can have a strong influ¬ 
ence. 

Upon including the mean held potential, the effects 
of interactions manifest within two distinct terms of 
the Boltzmann equation. The natural question arises 
whether or not there is some error akin to double count¬ 
ing due to the presence of both these terms. The collision 
term describes an instantaneous collision between exactly 
two particles within the gas, such that momenta is ex¬ 
changed between these two particles. This effect is en¬ 
tirely local, and occurs irrespective of the other particles 
in the gas. On the other hand, the mean held consists 
of a collective effect due to every single particle in the 
gas. In this sense the two terms are conceptually dis¬ 
tinct from one another. Serious problems begin to occur 
when the mean held interaction energy becomes partic¬ 
ularly signiheant (taking up a large fraction of the total 
energy in the gas). In such a situation, the collisions can 
begin to occur, not on the background of a translation- 
ally invariant potential energy landscape (as it is gener- 
ably assumed [14[) but rather on an appreciably varying 
potential energy landscape, caused by the mean held of 
nearby particles. These problems arise when typical val¬ 
ues of na 3 d approach or exceed unity. As we show below, 
this is not the case in our current realm of interest. 

In order to ascertain the relevance of the mean held 
in Eq. for our current simulation, we wish to con¬ 
sider the total mean-held energy per particle e m f in the 
gas, and compare this to the temperature. That is, we 
calculate 


where 


, . . 1 + 2x 2 3a7 2 arctanhVl — x 2 , . 

h{x) = -< 25 > 

is a function generally of order unity (although h{l) = 0 
since the angular average of Vdd is zero). In an attempt 
to draw some broad conclusions, we simply consider the 
prefactor in e m f and compare it to the temperature: 


1 N C dd 
kBT 48\/tt 3 o\oz ’ 


(26) 


In the experiment of Ref. fl5| which we are currently 
interested in, the quantity n is never more than k < 
0 .02, indicating that physics associated with the mean- 
held is likely to be insigniheant, at least to a hrst level of 
approximation. 

In other situations (involving higher densities, or larger 
dipole length scales), where k becomes appreciably large, 
incorporating the mean-held into the simulation may be 
necessary. The computational issues of doing so are, to a 
certain extent, manageable (see for example the vast lit¬ 
erature on particle-in-cell methods used to solve the ordi¬ 
nary Vlasov equation in the held of plasma physics [56j]). 
Briehy, the process involves binning the particles in po¬ 
sition space to hnd the density n(r, t), smoothing the 
density via convolution with a suitably chosen gaussian 
kernel, and then calculating t he p otential, using Eq. 
and ultimately the force F j56|. For issues relating to 
clarity, we currently wish to relegate further details of 
this procedure to a future publication. 


V. RESULTS FOR FERMIONS 


6mi = JN j rf3r n{r,t)U mf (r,t). 


( 22 ) 


We are only interested in placing an approximate upper- 
bound on the value of e m f, so we simplify the situation 
at hand by assuming the density of the gas (at any given 
time) is given by a gaussian distribution with cylindri¬ 
cal symmetry about the dipole-alignment direction which 
(solely for the purpose of this discussion) we assume to 
be along the 2 -axis; 


The choice of physical parameters in our simulation 
are taken directly from Ref. [l5j. These are 


N = 8 x 10 4 
T = 426nK 
m = 2.77 x 10“ 25 kg 
ad = 5.25nm 


i(r) = 


N 


(27t) 3 / 2 C 


■ exp 


± "z 





uj x = 2 tt x 393Hz 

x 2 + y 2 

z 2 1 

• (23) 

u4 0) = 2 t r x 38Hz > 

2<j\ 

2 cr 2 

lo z = 2n x 418Hz J 


total atom number 
initial temperature 
167 Er mass 
167 Er dipole length scale 

initial trap 


One could perform a more realistic calculation in the ab¬ 
sence of cylindrical symmetry, but analytic calculations 


Uamp — 14mS 
s = 1.8 


ramp time 
final trap, y —axis [see Eq. ©]■ 












7 


We vary the computational parameters Nt and AV un¬ 
til numerical convergence is achieved. This has typically 
occurred when Wp ft N, although we perform our sim¬ 
ulations right through to Nt = 4 x TV to thoroughly 
check the convergence. We find these simulations con¬ 
verge rather rapidly with AV [H3|, however we perform 
simulations right through to nAV = 0.35 (where n is the 
initial trap-averaged density), with Nt = 4 x TV, to be 
certain of convergence. 


A. Anisotropic pseudo-temperature 

To evaluate the rate of rethermalization, we find the 
standard-deviations of the test-particle distributions; for 
instance 


<7x{t) 



Nt 


2—1 


<7p* (*) 



N t 

Y^Pxi(t) 2 , 

2=1 


and equally for the y and z directions. We note that, a 
gaussian distribution provides a reasonably accurate ap¬ 
proximation to the instantaneous empirical distribution 
of test particles in the simulation. However, the moments 
above are well defined, regardless of whether this is the 
case or not. From these standard deviations, we can de¬ 
fine a time-dependent, anisotropic pseudo-temperature , 
related to the widths of the test-particle distribution 
function in each direction, relative to the instantaneous 
value of the trapping parameters, for instance; 


2 2 
_ mujx&x 

kn 


'T' _ Px 

Px mkB 


(27) 


and equally for the y and z axes. This definition makes 
particular sense in the case of a gaussian distribution. 
The two quantities; Tx and T Px above, can be combined 
into a single pseudo-temperature in the ^-direction (or 
in any direction) given by the mean; 


T x 



(28) 


The results of this analysis for the temperature along the 
z-axis is shown in Fig. [2j along with the experimental 
data of Ref. [li|- A more complete set of results, for the 
temperatures in all three directions is shown in Fig [3] 
An interesting observation we make is the apparent non¬ 
monotonic rethermalization behaviour of T z near p = 45° 
(this behaviour seems to exist right through 30° < /3 < 
60°). This behaviour was not observed in the experiment, 
likely due to the fact that it is a subtle effect which may 
be difficult to measure. Indeed we note in Fig. [5] (a), the 
scatter and error bars in the experimental data points 
appear to be of a similar size, or even larger than the 
magnitude of the non-monotonic hump in the theoretical 
result. 


B. Analyzing the rate-of-rethermalization as a 
function of p 


In order to define a rate-of-rethermalization it is cus¬ 
tomary to fit an exponential decay curve to the equilibra¬ 
tion dynamics shown in Figs.[2]and[3] For example, in the 
^-direction, one would write T z (t) = T z eq ^ + AT z e _t / Ti , 
where T z eo ^ (a fit parameter) is the equilibrated temper¬ 
ature, and T z eq ^ + A T z is the initial temperature (426nK 
in our case). The time-constant of this exponential decay 
curve, t z , is then written as 


Uz 

t z = — -r, 

no'F.B'y 


(29) 


where v = y/l6kBT / nm is the mean-velocity in the 
gas, and (Jf.b is the total cross-section of Eq. m av¬ 
eraged over all solid angles of the incoming relative mo¬ 
mentum p re i, such that dp = (327r/15 )a 2 d and ctb = 
8na 2 + (32-7r/45)a^. In this way, the quantity ndp^v rep¬ 
resents the mean-collision-frequency in the gas, and the 
quantity a can be conceptually thought of as the num¬ 
ber of collisions required for rethermalization. The exact 
same procedure can be applied to the x and y axes. In 
our current situation ot x , y , z will be a function the an¬ 
gle p between the dipole-alignment direction and the y 
axis. The results, which agree well with experimental 
data from Ref. [I5|, are shown in Fig. H) 

It should be noted that Refs. 3 151 ] compute a z in a 
simpler way, by approximating the short-time behaviour 
of the dynamics via the Enskog equation [58j]- This has 
also shown adequate agreement with the data, but gives 
considerably less detail than the present DSMC simula¬ 
tions. 


C. Trap-oscillations and covariances in position 
and momentum space 

The sudden change in the trap frequency along the y- 
axis gives rise to a breathing-mode along this direction 
(see Ref. [S9| for a discussion of this subject in the case 
of a classical gas with hard-sphere interactions). The 
oscillations are apparent in either the position variable 
T v , or the momentum variable l~ Py , but not in the sum 
T y which is plotted in Fig. [3] (since T Py and T y oscillate 
exactly out of phase with each other). This behaviour 
is shown in Fig. [5] The experiment of Ref. [3 nei¬ 
ther reported, nor searched for, any evidence of these 
oscillations or their damping periods (data was only an¬ 
alyzed along the z-axis). The frequency of the breathing 
mode is 2 uj y (t > tramp) mm- Collisions will eventu¬ 
ally cause this mode to damp out (intriguingly though, 
monopole modes are undamped in spherically symmetric 
harmonic traps). In order to quantify this, we subtract 
off the pseudo-temperature (shown by the red-dashed line 
in Fig. [5]), and fit a decaying sinusoid to the data; 

T y {t) — T y (t) « Ae -t//T080 sin [cot + 5] . (30) 















t (S) t (s) t (s) 


FIG. 2: (Color online) A comparison between the experimentally measured rethermalization process versus the results from 
the DSMC simulation. In each figure, the red solid line shows the result of the DSMC simulation, calculated analagously to 
Eq. (1281) . but along the z axis. The experimental data points are shown in blue with error bars. The agreement is reasonable, 
especially considering there was no post-processing made of the experimental data, nor any adjustments to the theory in order 
to produce these fits (no free parameters). 
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FIG. 3: (Color online) Shows the pseudo-temperatures along the x, y, and z axes (shown in red-dashed, blue-dot-dashed, and 
green-solid lines respectively) defined analogously to Eq. (1281) . In the first 14ms (the ramp time), the temperature along the 
y -axis increases as a result of the rapid change in the trap-frequency along this direction. After 14ms, the system relaxes down 
towards the new equilibrium state. The rate at which T x , T y and T z return to an equilibrium value displays a strong dependence 
on /3, which we explore further in section FV Bl An unexpected feature we observe is the non-monotonic path by which T z returns 
to equilibrium near ft = 45°. The effect is shown in greater detail in (1) by zooming in on the relevant part of (f). This is 
certainly an interesting consequence of the anisotropic dipole differential scattering, but note that the behaviour only occurs 
along one of the coordinate axes (the 2 -axis in this case) and, overall there is no violation of Boltzmann’s 17-theorem. 
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(5 (degrees) p (degrees) 

FIG. 4: (Color online) Shows a (the number of collisions re¬ 
quired for rethermalization) as a function of the angle (5 along 
the; (a) x (red dashed line) and y (blue solid line) directions, 
and (b) 2 direction. The data shown in (b) is taken from the 
experiment in Ref. ! 15] (data was not taken in the x and y 
directions). 


In the current experimental scenario the erbium gas lies 
firmly within the collisionless limit (trap frequency is sig¬ 
nificantly higher than the mean-collision frequency), and 
therefore the oscillation frequency is w = 2^1 + s Uy 0 ' 1 
i.e. twice the final trap frequency. Of course, if in¬ 
stead the experiment were in the hydrodynamic regime, 
rather than the collisionless regime, this would not be 
the case m m. We only fit to the region t > t ramp 
when the trap is no longer changing. The parameters A , 
T osc , and S are all fitting parameters. We then scale the 
time-constant r os c by the collision-frequency to give us 


C^OSC 

7"osc — — 

n(TF,BV 


(31) 


such that we can loosely interpret a osc as the number of 
collisions required for the breathing mode to damp out. 
Naively one might expect this to be the same as the a 
in section IV B1 and indeed we find distinct similarities, 
however the breathing mode takes considerably longer to 
damp out (a factor of 2 or more). The results for how 
a 0 sc depends on fd is shown in Fig. [G] note the qualitative 
similarity between Fig. 0 and the blue line in Fig. 0(a). 
We do not find that the other fitting parameters A and <5 
have any significant dependence on fd. However, A does 
depend on the size of the perturbation to the trap, and 
S depends on the ramp time t ramp (this is apparent in 
the instantaneous quench, for which analytic formulae 
are straight-forward). 

In contrast, breathing modes along the x and z axis 
are considerably less pronounced [wj. This is simply due 
to the fact that the perturbing force on the system in this 
situation is entirely along the y axis (see Fig. [[]). 

If the quench were performed instantaneously, a simple 
analytic solution is available in the extreme-collisionless 
limit: 


/(r, P, t) =/^ 2 b } [(x,z),{Px,Pz)\ x 


M exp 


~2(y Py 1 


y 

Py 


(32) 


where /mb' is the 2D Maxwell-Boltzmann distribution 
(along the x and z axes), A4 is a normalisation constant, 


and the covariance matrix 

*=($2) < 33 > 


is such that C = (j/ 2 ) - (j/) 2 , V = {yPy) ~ (y)(p v ), and 
8 = (p^) — (p y ) 2 . Note that f and 8 are proportional to 
the pseudo-temperatures T y and T Pv respectively, where 
as p is the covariance between position and momentum 
space. Ignoring collisions in the system, these variances 
evolve according to [61:] ; 


c = 

v = 
8 = 


Co 

2 L 

vm 


i + r + (i - r) cos (2 w^tj 


sin ^2 to 


r 1 / 2 _ r ~ 1/5 


8n 


i + r - 1 + (i-r- 1 )cos (2 w 


(34a) 

(34b) 

(34c) 


where Co = kBT/[m(ujy °' ) ) 2 ], and 8q = ksTm are the 
initial spatial and momentum variances (respectively), 

and F = (ujy^/ijjyA is the ratio of initial-to-final trap 


frequencies (squared). 

We have performed simulations of the cross¬ 
dimensional relaxation procedure in the case of an instan¬ 
taneous quench. The results are shown in Fig. 0 where 
we compare the simulation data to the analytic formulae 
of Eqs. (IM1) . The simulations reveal the increasing im¬ 
portance of collision-induced damping for times beyond 
several trap periods. The decay rate of the covariance p 
depends on the dipole angle fd. To within the numerical 
accuracy of these simulations, we find that the rate at 
which p decays, and the dependence this decay has on 
fd, is extremely close to that for <f and 8 (the pseudo¬ 
temperatures) shown in Fig. [6] 


D. Quantum many-body effects 

The Boltzmann equation, as written in Eq. m , treats 
the many-body dynamics of the system entirely in terms 
of classical mechanics. For our comparison with the ex¬ 
periment in Ref. [l5j , this may conceivably be a source of 
error. In 1928, Nordheim made adjustments to the Boltz¬ 
mann equation to account for the quantum-mechanical 
effects of Fermi-blocking and Bose-enhancement 62|. 
The net result of Nordheim’s work was an alteration to 
the collision integral: 

Cn[/] = /^/^S IP - Pi|X 

[ff[ (1 ± h 3 f) (1 ± h 3 h)-ff 1 (1 ± h 3 f) (1 ± h 3 f[)] 

(35) 

where h is Planck’s constant, and the + sign applies to 
identical bosons (Bose enhancement) while the — sign 
applies to identical fermions (Fermi blocking). From 
this point of view, the quantum many-body effects in 
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FIG. 5: (Color online) The rapid change in trapping frequency 
along the y -axis generates a large breathing mode along this 
direction. This is shown above in plots of T y versus time for a 
variety of different values of p. These breathing modes exist 
also in the momentum distribution, T Py , and look identical 
to the plots above except that the oscillations are exactly 7r- 
radians out of phase (leading to the monotonic behaviour in 
T y shown in Fig. dj. The dashed (red) line in each of the 
figures is T y . We use Eq. J30J as a fit to the decay of this 
breathing mode. The breathing mode dynamics along the x 
and z axes are barely noticeable in our simulations. 



FIG. 6: The breathing mode along the y -axis is damped over 
a time-scale r os c found from Eqs. (1301) and (1311) . The depen¬ 
dence on p is shown above. Note the qualitative similarity 
a OS c (shown above) has to a y in the blue line of Fig. [4] (a). 
However, the oscillations take considerably longer to damp 
than the envelope, as a osc > a s . 


the system are determined by the phase-space density 
(see Ref. [63] for a discussion, and recent results, on the 
fermionic gas). Specifically how many particles occupy 
a volume of phase space equal to h 3 . If this number is 
much less than one, quantum effects should be small, if 
this number is comparable to one, quantum effects will 
be important. The maximum phase-space density is plot¬ 
ted in Fig. [8] as a function of time for two different values 
of p. From this, we conclude that quantum many-body 
effects will have a negligible effect on the dynamics at 



FIG. 7: (Color online) Comparison between the DSMC simu¬ 
lation for an instantaneous quench and the analytic formulae 
in Eqs. ([34]). In (a) we plot the covariance between position 
and momentum space, and in (b) we plot the variance in po¬ 
sition space (proportional to Ty). This particular data is for 
a dipole alignment direction of ft = 0. The DSMC simula¬ 
tion is shown by the solid (blue) line, the analytic formulae 
by the dashed (red) line. The analytic formulae do an excel¬ 
lent job of correctly predicting the amplitude and phase of 
the oscillations. For this ratio of collision-to-trap frequency, 
the damping becomes appreciable on the order of several trap 
periods. 




FIG. 8: Shows the (maximum) number of particles in a vol¬ 
ume element of phase space equal to h 3 as a function of time 
for two separate dipole-alignment angles; (a) ft = 0, and (b) 
P = 45. The phase space density decreases as the system equi¬ 
librates to a higher final temperature. From this, we estimate 
that quantum many-body effects are indeed small enough to 
be neglected (at least as a first approximation). 

this temperature. This goes some way in explaining the 
reasonably good agreement between our theory and ex¬ 
periment in this case. We do not expect our theory to 
provide quantitative accuracy at significantly lower tem¬ 
peratures, although modifying our algorithm to account 
for the mechanism of Bose-enhancement/Fermi-blocking 
is a future goal of this project. Speculating further on 
this, we note that the Boltzmann-Nordheim equation will 
have, not only a (potentially) different path to equilib¬ 
rium, but also (at lower temperatures) a different equilib¬ 
rium state as well (the famous Bose-Einstein and Fermi- 
Dirac distributions). How this would affect the depen¬ 
dence of ot x ,y,z on P is an interesting and open question. 

VI. RESULTS FOR BOSONS 

It is very straightforward to repeat these simulations 
for a system of bosons simply by replacing gp with in 
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the differential scattering cross-section and ctf —► ob (see 
Eqs. (Illal) . (Illbl) . and m in section [Til Cl) . We choose 
to keep the geometry of the trap, the atomic species, 
and the number of particles the same as that which was 
used in section m for fermions. We set the s-wave scat¬ 
tering length a = 0, to emphasize the peculiarities of 
the anisotropic dipolar differential scattering. The dis¬ 
tinctions between bosonic versus fermionic scattering be¬ 
haviour naturally alters details of the rejection sampling 
algorithm (see appendix 0 and changes the results, but 
there is no conceptual change in what we are doing, so we 
provide less detail than we did for fermions. In addition, 
experimental data does not yet exist for bosons, so we 
cannot make the same comparisons in that respect. 

Figure [9] shows the rethermalization of the pseudo¬ 
temperatures for bosons (analagous to Fig. [3] for 
fermions). Somewhat ironically, in the context of low- 
energy scattering, the rethermalization procedure takes 
approximately three times longer for bosons than for 
fermions with the same density and dipole-moment. This 
is due to the factor of three difference (for a = 0) be¬ 
tween the angularly averaged total cross sections ctf and 
ctb 0- Increasing the s-wave scattering length a would 
naturally change this situation. The nature of the differ¬ 
ential cross-sections are such that a nonmonotonic rether¬ 
malization process is not observed for bosons [as it was 
in Fig. [3] (1)]. Figure [TUI fal-fcl shows the number of col¬ 
lisions required for rethermalization as a function of /3. 
In (d) we show the maximum phase-space density as a 
function of time for the case /3 = 30°. Again this in¬ 
dicates that the Boltzmann equation should provide an 
approximately accurate theoretical description at these 
densities and temperatures. Figure [11] shows the num¬ 
ber of collisions required to damp out the breathing mode. 
Note the qualitative similarity between a osc h 1 Fig. fill 
and a y in Fig. [TO] (b), but with a quantitative difference 
of approximately a factor of two. 


VII. CONCLUSIONS AND DISCUSSION 

In this article we have developed a DSMC numerical 
algorithm to solve the Boltzmann equation for an ultra¬ 
cold dipolar gas. We have used this method to study the 
cross-dimensional relaxation dynamics of a dipolar gas 
via a full simulation of the phase-space dynamics. Where 
applicable, we have compared our numerical results with 
the experimental data of Ref. |l5( and found favourable 
agreement. This suggests that the DSMC algorithm pro¬ 
vides a quantitative method for understanding the nor¬ 
mal component in a dipolar gas. This is a promising 
result. The method is suitable for both fermions and 
bosons, although experimental data currently exists only 
for fermions. The method and results direct our inter¬ 
est toward a new regime where interactions in the gas 
manifest from collisions rather than the mean-field. 

More specifically, we have studied the damping of trap 
breathing modes in the system and quantified the pro- 



t (s) t (s) 

FIG. 9: (Color online) Shows the pseudo-temperatures along 
the x, y, and z axes (shown by red-dashed, blue-dot-dashed, 
and green-solid lines respectively) as a function of time for the 
bosonic dipole scattering cross-section. Experimental data 
has not been taken for this case, but we observe the rethermal¬ 
ization rates showing a strong dependence on /? (particularly 
along the x-axis). 




135 180 



FIG. 10: (a), (b), and (c) show a x , a y , and a z (the number 
of collisions required for rethermalization) as function of j3 
in the case of bosons, (d) shows the maximum phase space 
density as a function of time for /3 = 30°. 



FIG. 11: The decay of the breathing mode along the j/-axis 
in the case of bosons. Again, there is a strong qualitative 
similarity between this curve and the curve in Fig. [TO] (b), 
but an important quantitative difference in that aw is larger 
by approximately a factor of two. 
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nounced dependence of rethermalization on the dipole- 
alignment direction. We find the breathing mode takes 
significantly longer (approximately a factor of two) to 
decay than the envelope for rethermalization, which is 
found by averaging over momentum-space and real-space 
dynamics. 

Our current work is entirely focused on the thermal 
gas, above quantum degeneracy. There are several rea¬ 
sons why understanding this normal component of an 
ultra-cold dipolar gas is important. For instance, at¬ 
tractive interactions along the dipole alignment direc¬ 
tion (due to the mean-field) can destabilise the sys¬ 
tem [ 2 , iSJei. Thermal energy can counter-act this 
instability [5j, [ 66 ], therefore we expect the normal com¬ 
ponent to have a qualitative, as well as quantitative, role 
in the dynamics. Our method presented here, if com¬ 
bined/coupled with one of the many low-temperature 
theories (e.g. [13, 68 |) would constitute a complete fi¬ 
nite temperature description of dipolar gases (in the same 
vein as the Zaremba-Nikuni-Griffin formalism of regular 
Bose-condensates [43J, |&9, 7CJ]). This remains as work-in¬ 
progress. 

The method used in this paper (DSMC) is a remark¬ 
ably versatile tool, potentially capable of simulating a 
multitude of out-of-equilibrium scenarios. Extending it 
into a regime where many-body quantum mechanical be¬ 
haviour becomes prevalent (beyond the simple two-body 
scattering level which plays such a vital role in our cur¬ 
rent work) is a direction which we intend to take this 
research. Possible avenues for doing so include, incorpo¬ 
rating the effects of Bose-stimulation and Pauli-blocking 
into the differential scattering cross sections, as pre¬ 
scribed by Nordheim (Hj, see Eq. (1551) . This requires 
modifications to the DSMC algorithm, which were orig¬ 
inally introduced in the context of nuclear equations of 
state, particularly during heavy ion collisions EM Q- 
The basic ideas have seen application in ultra-cold atomic 
systems of fermions, see Refs. EzHil- Another possibil¬ 
ity, perhaps more relevant for bosonic systems, involves 
coupling the Boltzmann equation (the purely classical 
version may suffice) to an equation describing the super¬ 
fluid component in the system. For example one could 
consider using the well-known Gross-Pitaevskii equa¬ 
tion [[70], or the more sophisticated c-field techniques [67f] ■ 
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Appendix A: Rejection sampling algorithm 

The procedure of rejection sampling is not new Jdj, 
but for completeness, we provide a brief description of 
the details specific to our situation. A more thorough 
description of the algorithm in general can be found in 

Ref. (zl. 


1. Fermions 

To sample from Pp{8 , <jr, rj) defined in Eq. (fl3l) . the 
strategy is to start from a simpler distribution (which 
is easy to sample), call it g(9,cj)) = l/(27r 2 ), and (appro¬ 
priately) reject those samples which were unlikely (recall 
that we only need to sample 9 and 4> since rj is given to 
us by the (already known) incoming relative momentum 
of the collision pair). The algorithm goes as follows: 

1. Sample (9,(f>) from g(0, <(>), and sample u from 
7/(0,1) (the uniform distribution over the unit in¬ 
terval) . 

2. Check whether u < Pp{9, (j)\ rj)/[Mg(9, ^)] 

where M is an upper-bound such that 

M > Pp(9 , </>; rj)/g{0, <f> ) for all 9 and <f>. 

3. If step 2 holds true, accept ( 9 , q i) as a realisation of 
Pp. If it does not hold true, reject (6, (j>), and begin 
over at step 1 . 

In order to find the upper bound M(rj) we transform 
to the collision-reference-frame, where 


Pf(9,^>-,v) 


6 sin( 0 ) [cos( 0 ) (cos 2 (rj) — cos 2 ((/>) sin 2 ( 77 )) + cos(^) sin( 0 ) sin( 2 ry )] 2 
7 T (3 + 18 cos 2 (rj) — 13cos 4 (?y)) 


(Al) 


Using standard optimisation methods, we find the max¬ 
imum value of Pf(9, <j>\ rj) occurs at (f > max = 0, and 
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^max — < 


( \/7+cos(4r ) )- x /2 sin 2 (T))(17-cos(4^)) \ . /o ^ ^ Q /A 

acos | -- 9^/3 - I 9 < 7r /4 or v 2 < V < 37T/4 

7+cos(4r;) —-\/2 sin 2 (r;)(17 —cos(4t;)) 


(A2) 


acos — 


2\/3 


7 r /4 < 77 < 7 r /2 or 37 r /4 < 77 . 


r 


From this, we define M = 2t: 2 P^ na ^\rj) where, 
p(max)/ \ _ 6 cos ($ max 2 ?y) sin($ max ) 


2. Bosons 


j (n) = """" rmal _ ““*v " max 7 t A 3 l The procedure for bosons is essentially equivalent, ex- 

w tt( 3+ 18cos 1 2 3 77- 13cos 4 5 6 7 8 9 10 11 r?)‘ 1 J cept with; 




2 sin(0) [—2 + 3 cos 2 ( 77 ) + 3 cos 2 (</> ma x) sin 2 (ry)] ‘ 
7 r (11 — 30 cos 2 (t 7) + 27cos 4 (^)) 

I- 


(A4) 


^max — 7 t/2, 


0 77 < atan(-\/ 2 ) or 77 > 7 r — atan(v / 2 ) 

7t/ 2 atan(-\/2) < 77 < 7 r — atan(-\/2) . 

(A5) 


and 

nimaxj/ \ 

P B '(7/) = 


(max) ( 2 [-2 + 3 cos 2 r] + 3 cos 2 0 max sin 2 77 ]" 


7 r (11 — 30 cos 2 ?7 + 27 cos 4 rj) 


(A 6 ) 


Note that P B in Eq. (IA4I) factorizes into a product of two 
functions involving only 6 and only cj>. This was not the 
case for the fermionic cross section, see Eq. ED. This 
allows for the sampling algorithm to be more efficient in 
the case of bosons than it is for fermions, since 9 can be 
sampled directly. 
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